The itinerant ferromagnetic phase of the Hubbard model 
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Using a newly developed quantum Monte Carlo technique, we provide strong evidence for the sta- 
bility of a saturated ferromagnetic phase in the high-density regime of the two-dimensional infinite-t/ 
Hubbard model. By decreasing the electron density, a discontinuous transition to a paramagnetic 
phase is observed, accompanied by a divergence of the susceptibility on the paramagnetic side. 
This behavior, resulting from a high degeneracy among different spin sectors, is consistent with an 
infinite-order phase transition. The remarkable stability of itinerant ferromagnetism renews the hope 
to describe this phenomenon within a purely kinetic mechanism and will facilitate the validation of 
experimental quantum simulators with cold atoms loaded in optical lattices. 
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Ever since classical antiquity, ferromagnetism has at- 
tracted the attention of natural philosopers. [l| A proper 
understanding of this phenomenon was only made pos- 
sible by the advent of quantum mechanics, from the 
early interpretations [2|, |3| to its modern realizations in 
quantum simulators engineered by means of cold atomic 
gases, y In some solids, such as transitions metals, the 
spin-independent nature of interactions has led to conjec- 
ture that long-range magnetic order is due to an itinerant 
mechanism in which the Coulomb interaction and the 
Pauli exclusion principle play a fundamental role. The 
single-band Hubbard model, possibly the simplest and 
most studied lattice model of correlated electrons, was 
first thought to encompass a minimal description of itin- 
erant ferromagnetism. [5| Recent experiments on ultra- 
cold atoms have shown that a gas of fcrmions may become 
ferromagnetic because of repulsive interactions. [4 This 
important result and subsequent numerical calculations 
in the continuum [6| suggested that this phenomenon is 
universal and independent upon the details of the inter- 
action, thus renewing the interest in the Hubbard model 
with a large Coulomb repulsion, U . In spite of its simplic- 
ity, exact solutions of the Hubbard model are not avail- 
able in more than one spatial dimension, leaving the ques- 
tion of the stability of a ferromagnetic phase unsolved. 
One of the very few exact results that is known is due to 
Nagaoka, [Tf who proved a theorem stating that, in the 
infinite- C/ limit, a single hole stabilizes a fully polarized 
ground state. Following this pioneering work, much effort 
has been devoted to s tudy the fully-polarized state for fi- 
nite hole densities. |8l4l4|- However, a general consensus 
on the stability of ferromagnetic phases is still lacking. 

In this Letter we present new results for the infinite-C/ 
Hubbard model, based on accurate fcrmionic quantum 
Monte Carlo (QMC) simulations, which indicate that at 
high electron density the Nagaoka state is stable not only 
with respect to the paramagnetic phase, but also with 
respect to other previously proposed partially polarized 
states. [12] A non-trivial transition to a paramagnetic 
phase is observed upon decreasing the electron density. 



Near the transition this phase is characterized by highly 
degenerate states with different values of the total spin, 
thus indicating a divergence of the magnetic susceptibil- 
ity, consistent with an infinite-order phase transition. [l5| 
The QMC simulation of systems of interacting elec- 
trons is beset by the antisymmetry of the ground-state 
wave-function which, at variance with bosons, prevents 
it from being treated as the stationary distribution of a 
diffusion process. The main attempt to cope with the en- 
suing difficulties is the so-called fixed-node (FN) approx- 
imation, which, for lattice models, amounts to defining 
an effective Hamiltonian whose ground state energy is a 
variational upper bound to the exact energy, [iq If com- 
plemented by an accurate variational ansatz for the wave- 
function, the FN method provides a method to study the 
properties of large fermionic systems making possible re- 
liable extrapolations to the thermodynamic limit. Unfor- 
tunately, the nature of the approximation does not allow 
for an estimate of the residual error, which not rarely can 
lead to biased results. However, the infinite-C/ Hubbard 
model belongs to an interesting class of Hamiltonians 
whose eigenstates of fermionic symmetry are sufficiently 
close in energy to the bosonic ground state, to allow 
them to be treated on an equal footing; for this class of 
Hamiltonians we propose a strategy to overcome the sign 
problem via the dissection of the excitation spectrum of 
the corresponding bosonic auxiliary problem, providing 
an essentially unbiased scheme for medium-size fcrmionic 
systems. 

Fermionic- correlation method. The spectrum of a 
Hamiltonian of identical particles, "H, can be classified 
according to the irreducible representations of the sym- 
metric (permutation) group. The Pauli principle states 
that only totally antisymmetric states are physically al- 
lowed for fcrmions, but mathematical states of any sym- 
metry can also be considered. In particular, the (un- 
physical) state of lowest energy is in general totally sym- 
metric, so that the fermionic ground state can be for- 
mally considered as an excited state of a bosonic system. 
As such, it can be studied via excited-state techniques, 



provided the Bosc-Fcrmi gap is not too large with re- 
spect to the physical gap in the fermionic sector of the 
spectrum. Let j^P^) be the bosonic ground state of the 
system and A an arbitrary observable. A recent exten- 
sion of the reptation QMC method [17| to lattice mod- 
els [l8[ allows for an efficient and unbiased evaluation 
of imaginary-time t = it correlation functions, Cy({T) = 
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the Heisenberg representation of A. 

The connection of such correlation functions with the 
excited states 1^*^) of H is obtained by considering the 
Lehman spectral representation, 
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where Ak = E^ — E^ are excitation energies with respect 
to the bosonic ground state. Selection rules act in such a 
way as to exclude from Eq. ([1]) those excited states whose 
symmetry is different from that of the state ^|\l/°). In 
particular, if A is chosen to be totally antisymmetric with 
respect to permutations, only fermionic (ground and ex- 
cited) states would contribute to C_A.iT). For example, if 
A is the local operator whose coordinate representation 
is the ratio between the fermionic and bosonic ground- 
state wave-functions {Af{n) = (n|\I'j)/(n|$°), where \n) 
denotes the many-body lattice configuration), the corre- 
lation function Cj[{t) would be proportional to the single 
exponential e"^""^. 

In practice, neither the bosonic nor the fermionic 
ground state are known exactly and only variational ap- 
proximations to them are available, which we denote 
here by |$b) and |<&/), respectively. Correspondingly, 
the antisymmetric observable is defined as Af{n) = 
{n\^ f) / {n\^i,) . In this way, the leading coefficient of the 
expansion is given by (\E'°|.A/|\[''^) ~ ($/|^^) and can be 
systematically maximized improving the quality of the 
variational states. The energy of the fermionic ground 
state can be then extracted either directly by noticing 
that E^f = E^ — limT-_>oo [9r logC^(r)] or, indirectly, by 
fitting the exponential decay of the correlation function 
of Eq. dl]) and extracting the smallest energy gap. In or- 
der for this procedure to make any sense, it is necessary 
that the (unphysical) Bose-Fermi gap is not too large 
with respect to the physical excitation energies in the 
fermionic sector of the spectrum. If this condition is not 
met, the anti-symmetric correlation function gets effec- 
tively extinguished before the selection of the fermionic 
ground state from its excitation background is attained 
by imaginary-time evolution. This condition is actually 
verified for infinite-C/ fermionic Hubbard models of mod- 
erate size, where the properties of the system are not too 
dissimilar from those of a system of hard-core bosons. 
The condition of a small fermion-boson gap is also met 
in other interesting systems, where the effects of statis- 
tics on the total energy are overwhelmed by the effects 
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Figure 1: Energy E'J{t) = E^ — 9^ log C^(r) as a function 
of the imaginary-time r for L — 50 and N = 42 electrons 
and different magnetizations. The dashed horizontal lines are 
FN energies, while the solid lines are the energies as obtained 
fitting the imaginary-time correlations. 



of correlations, such as the low-density electron gas, liq- 
uid '^He, quasi-unidimensional systems and mixtures of 
bosons and fermions. 

We remark that even in the most favorable cases all 
these considerations can only hold for systems of moder- 
ate size because the Fermi-Bose gap is an extensive prop- 
erty, whereas the physical gap in the Fermi sector of the 
spectrum is intensive, being determined by quasi-particle 
effects. This fermionic-correlation method is related to 
the transient estimate (TE) method for the fermionic 
ground state, [19| or its generalization for a few excita- 
tions. [20[ However, TE works with ratios of decaying cor- 
relation functions, thereby reducing the signal/noise ra- 
tio, and typically uses sub-optimal bosonic guiding func- 
tions, with increased fluctuations in the weights of the 
random walks. 

The calculation of spectra from imaginary-time cor- 
relation functions is in general an ill-posed problem. 
In practice, however, sharp peaks with strong spectral 
weight can be reliably extracted if the correlation func- 
tion is known with good statistical precision for suffi- 
ciently long times. [17|,|21| In the present work this con- 
dition is met even for systems of several tens particles, 
due to the relatively small gap between the fermionic and 
bosonic ground states, as well as to the good quality of 
the variational states. 

The model. The Hamiltonian of the infinite-f/ Hub- 
bard model reads: 



-Hf = -t 



E 



VgcI^c.^Vg + h.. 



(2) 



where cj ^ (ci.a) creates (destroys) an electron on site i 
with spin a; {i,j) denotes nearest-neighbor site pairs and 
the Gutzwiller projector Vg forbids double-occupancy. 
In the following, we will consider a square lattice and 
take t = 1 as the energy scale. The total number of sites 
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Figure 2: Ground-state magnetization of the infinite-f/ Hub- 
bard model on tiie square lattice. The shaded area represents 
a small region of uncertain attribution due to the effect of the 
residual Monte Carlo error. 



will be denoted by L and the number of electrons by 
A'^. In the following, we present the results for different 
magnetizations m = {n^ — n^)/ [n^ + n^) and densities 
n = n^ +ni. 

Relatively sirnple variational wave functions have been 
constructed, [9|,ll0| by flipping one (say up) spin with re- 
spect to the saturated ferromagnetic state. The flip of the 
spin leads to a gain of kinetic energy for the down spin, 
but also a loss in the spin-up kinetic energy (since the 
motion of spin up electrons is restricted by the necessity 
of avoiding double occupancy). Here, we consider 

{n\<Pf) =. Jf{n) X Det [MR])} x Det [MR])} , (3) 

where the Jastrow factor J^[n) = exp X^i i ^ii"'i"j 
multiplies two Slater determinants that are constructed 
by applying backflow correlations to single-particle or- 
bitals for up and down spins. [22| The correlated orbitals 
are defined by (/)fc(i?J) = 0^(i?J) + 6^ Eh,,.' '^"(^z"'), 
where bk are orbital-dependent backflow parameters, 
0°(i?J) are plane waves, and the sum includes all nearest 
neighbors of the j-th particle, thus preserving the spin 
rotational invariance. The proposed backflow wave func- 
tion (j3]) encodes the effect of correlation on the deforma- 
tion of the free-orbitals nodal structure and consistently 
catches much of the physics of previous treatments, [9|- 
ll| while leaving room for a systematic improvement with 
the QMC methods. 

The bosonic counterpart of the model studied is a 
purely kinetic hard-core bosons Hamiltonian, where the 
fermionic operators are substituted by bosonic ones. Our 
QMC method is particularly suitable to study the high- 
density region, namely few holes close to full filling {N = 
L), where the boson- fermion gap is very small and in- 
creases upon decreasing the density. |23| The bosonic trial 
state is given by a Jastrow wave function (n|$b) = j''{n), 
which is similar to the fermionic one (but with different 
parameters T^ „) and represents an excellent ansatz for 
the bosonic ground state. [2J] In all cases, the variational 
parameters are fully optimized minimizing the variational 
energy with the method of Ref. |25 . 

Results. The fermionic correlation technique remains 



efficient up to relatively large system sizes (i.e., L ~ 
50 -^ 100) and allows us to reach numerical results, which 
are exact within statistical accuracy. In Fig. [TJ we re- 
port our results for L = 50 and N = 42 electrons, for 
different values of the magnetization, m. In addition, we 
also report the results based upon the FN approach. The 
possibility to obtain numerically exact results on rather 
large systems allows us to assess the accuracy of the FN 
method that can be extended to much larger sizes (i.e., 
L < 1000), without any numerical instability. Thanks 
to backfiow correlations, we get a considerable improve- 
ment u pon the standard plane waves that has been used 
in Ref. 
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There is a small difference between the FN 
results and the energies obtained by the imaginary-time 
correlations, indicating a very small residual FN error, 
namely AE/t < 0.01. 

In Fig. [21 we report the overall phase diagram obtained 
by considering large-scale FN calculations. A saturated 
ferromagnetic phase is stable for n > 0.75, while for 
smaller densities a paramagnetic ground state is found. 
The narrow shaded region denotes the incertitude due 
to the residual numerical error, which can be estimated 
by comparing the FN energies with the exact ones (ob- 
tained from the fermionic correlations) on smaller clus- 
ters, see Fig.[TJ This direct comparison puts us on secure 
grounds as concerns the robustness of the dependence of 
the ground-state magnetization on the electron density. 

In Fig. 13] we display the dependence of the ground- 
state energy upon magnetization, for different values of 
the electron density. The remarkable feature emerging 
from this figure is the strong fiattening of the energy as a 
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Figure 3: Difference between the energy per site of different 
magnetizations and the one of the saturated ferromagnet as a 
function of the density n. The cases with L = 200 (squares) 
and L = 400 (circles) are reported; lines connecting points 
are a guide for the eye. 



function of the magnetization (i.e., the spin) close to the 
transition between the fuUy polarized ferromagnet and 
the paramagnetic state. Indeed, at low and high densi- 
ties the energy has a monotonic behavior as a function 
of the magnetization m. At low density a clear minimum 
exists at TO = 0, typical of a paramagnetic phase, where 
the curvature of the energy-versus-magnetization curve 
witnesses to a finite spin susceptibility. On the other 
hand, in the high-density ferromagnetic phase, E{m) dis- 
plays a well defined minimum for m = 1. By approach- 
ing the transition, E{m) becomes flatter and flatter, sug- 
gesting that the susceptibility may diverge at the critical 
point. Although we cannot exclude a tiny region with 
a finite but non-saturated magnetization, these results 
would suggest that the paramagnetic-to-ferromagnetic 
transition is not due to a simple level crossing, namely to 
the creation of a local minimum in E{m) at to = 1 that 
eventually prevails over the paramagnetic one, but rather 
to the progressive flattening of the whole E[m) curve. 

Our scenario is compatible with an infinite-order phase 
transition, which, in general, is described by E{m) = 
(g — gc)m'^ + bm?^ , where r — > oo; a phase transition is 
obtained by varying the order parameter g (in our case 
the electron density) across its critical value gc- The crit- 
ical exponent of the order parameter is /3 = l/(2r — 2), 
generating a jump from zero to the saturation value for 
r ^ oo. Moreover, the susceptibility x ^ ^±l\g ^ .9c T 
has an exponent 7 = 1 independent of r, with an ampli- 
tude ratio A_/A_|_ that vanishes for r — >■ 00. [15[ Even 
though the order parameter shows a finite jump, like in 
ordinary first-order phase transitions, there is no hys- 
teresis. We have indeed verified that the ground-state 
energy is a convex function of the electron density, im- 
plying a finite compressibility in the neighborhood of the 
ferromagnetic-paramagnetic transition. This picture im- 
plies that spin-fiip excitations over the fully polarized 
state arc non-interacting at the transition point. In fact, 
we find that, at small distances, the minority spins repel 
each other, whereas at large distances they do not inter- 
act. In the variational wave function, this fact generates 
a sizable repulsive short-range Jastrow factor, while at 
long range the V^ ■ pseiidopotential vanishes. 

Conclusions. In this Letter we have analyzed with high 
accuracy the magnetic phase diagram of the fcrmionic 
Hubbard model on the square lattice in the limit of in- 
finite on-site repulsion U. By the combination of differ- 
ent QMC methods, we are able to give a very precise 
determination of the transition between the ferromag- 
netic and the paramagnetic states. Interestingly, all spin 
excitations become essentially gapless at the transition, 
possibly indicating that the transition is of infinite or- 
der. [l5| Compared to previous calculations, this is the 
first time that such a behavior is observed. Indeed, given 
the extreme difficulty to treat this highly-correlated sys- 
tem, most of the theoretical efforts were limited to study 
very high densities or a single spin fiip. [8J-UJ| 



Our results pave the way to a better understanding 
of itinerant ferromagnetic phenomena in both traditional 
condensed matter systems and in recent and forthcoming 
realizations of such phases in cold atomic gases. Indeed, 
the recent achievements for the realization of interact- 



ing fermionic systems trapped in optical lattices [26[ will 
most likely lead to experimentally probe the strongly- 
correlated regime of the Hubbard model at sufficiently 
low temperatures. Finally, the generality of the numeri- 
cal methods introduced in this Letter will also offer new 
insights in other strongly correlated fermionic systems 
where currently available analytical and numerical treat- 
ments may fail to offer a quantitative or even qualitative 
account of the relevant physical properties. 
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